{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018  by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time-domain viscoelastic equations of motion \n",
    "\n",
    "In the last lesson we proved that the Generalized Maxwell model (GMB) is able to describe a constant Q-spectrum, necessary for realistic SH wave propagation in viscoelastic media. Now, we only have to transform the frequency domain GMB stress-strain relation back to time-domain and assemble the forward problem in a form which can be solved by finite-differences.\n",
    "\n",
    "## GMB stress-strain relation: from frequency to time domain\n",
    "\n",
    "By combining a Maxwell Standard Linear Solid (SLS) in parallel with additional Maxwell bodies, we derived the **Generalized Maxwell Body (GMB)**\n",
    "\n",
    "<img src=\"images/GMB.png\" width=\"50%\">\n",
    "\n",
    "which can be described in the frequency domain by the stress-strain relation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\tilde{\\sigma} = \\tilde{\\mu}_{GMB} \\tilde{\\epsilon} \\notag\n",
    "\\end{equation}\n",
    "\n",
    "with the complex shear modulus:\n",
    "\n",
    "\\begin{equation}\n",
    "\\tilde{\\mu}_{GMB} = \\mu_0 + \\delta \\mu \\sum_{l=1}^{L} \\frac{i a_l \\omega}{i \\omega + \\omega_l}. \\notag\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mu_0$ denotes the relaxed shear modulus, $\\delta \\mu$ the modulus defect, L the number of Maxwell bodies, $a_l, \\omega_l$ weighting coefficients and relaxation frequencies of the l-th Maxwell body to achieve a constant Q-spectrum, while $\\omega$ is  the circular frequency within the frequency range of the source wavelet.\n",
    "\n",
    "For time-domain FD modelling, we obviously need to transform the complex modulus to time-domain by inverse Fourier transform leading to the relaxation function $\\Psi(t)$\n",
    "\n",
    "\\begin{equation}\n",
    "\\Psi(t) = \\frac{1}{2\\pi}\\int_{-\\infty}^{\\infty}\\tilde{\\mu}_{GMB} \\exp(i\\omega t) d\\omega. \\notag\n",
    "\\end{equation}\n",
    "\n",
    "This leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\Psi(t) = \\biggl\\{\\mu_0 + \\delta \\mu \\sum_{l=1}^L a_l \\exp(-\\omega_l t) \\biggr\\} \\cdot H(t), \\notag\n",
    "\\end{equation}\n",
    "\n",
    "where $H(t)$ is the Heaviside step function. Using the definition of the modulus defect as the difference between the unrelaxed $\\mu_u$ and relaxed shear modulus $\\mu_0$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\delta\\mu = \\mu_u - \\mu_0\\notag\n",
    "\\end{equation}\n",
    "\n",
    "we can replace the relaxed by the unrelaxed modulus ([Emmerich & Korn 1987](https://library.seg.org/doi/10.1190/1.1442386), [Moczo & Kristek 2004](http://www.spice-rtn.org/library/publications/publicationbook.2007-05-08.3817028446.html)):\n",
    "\n",
    "\\begin{equation}\n",
    "\\Psi(t) = \\biggl\\{\\mu_u - \\delta \\mu \\sum_{l=1}^L a_l\\biggl(1 - \\exp(-\\omega_l t)\\biggr) \\biggr\\} \\cdot H(t). \\notag\n",
    "\\end{equation}\n",
    "\n",
    "Before inserting the relaxation function into the viscoelastic stress-strain relation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma = \\int_{-\\infty}^{t}\\dot{\\Psi}(t-t')\\epsilon(t')dt' = \\dot{\\Psi} * \\epsilon\\notag\n",
    "\\end{equation}\n",
    "\n",
    "we have to take its first time derivative:\n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{\\Psi}(t) = - \\delta \\mu \\sum_{l=1}^L a_l \\omega_l \\exp(-\\omega_l t) \\cdot H(t) + \\biggl\\{\\mu_u - \\delta \\mu \\sum_{l=1}^L a_l\\biggl(1 - \\exp(-\\omega_l t)\\biggr) \\biggr\\} \\cdot \\delta(t), \\notag\n",
    "\\end{equation}\n",
    "\n",
    "with the Delta function $\\delta(t)$, where we used\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial H(t)}{\\partial t} = \\delta(t)\\notag\n",
    "\\end{equation}\n",
    "\n",
    "Inserting the time-derivative of the relaxation function $\\dot{\\Psi}$ into the viscoelastic stress-strain relation yields:\n",
    "\n",
    "\\begin{align}\n",
    "\\sigma(t)=& -\\int_{-\\infty}^{t}\\delta\\mu \\sum_{l=1}^L a_l \\omega_l \\exp(-\\omega_l (t-t'))\\cdot H(t-t')\\; \\epsilon(t') dt'\\notag\\\\\n",
    "&+ \\int_{-\\infty}^{t} \\mu_0 \\cdot \\delta (t-t') \\cdot \\epsilon(t') dt'\\notag\\\\\n",
    "&-\\int_{-\\infty}^{t}\\delta\\mu \\sum_{l=1}^L a_l \\biggl(1 - \\exp(-\\omega_l (t-t'))\\biggr)\\cdot \\delta(t-t')\\; \\epsilon(t') dt'\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "Evaluating the Heaviside and Delta functions leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma(t) = \\mu_u \\cdot \\epsilon(t) - \\delta\\mu \\sum_{l=1}^L a_l \\omega_l \\int_{-\\infty}^t\\epsilon(t')\\cdot \\exp(-\\omega_l (t-t')) dt'\\notag\n",
    "\\end{equation}\n",
    "\n",
    "The problem with this expression is the remaining time-convolution, which would be required to apply at each spatial grid point and at each time step during the FD time stepping, resulting in significant computational costs. To get rid of this convolution, we use a standard physicist trick: Step 1, we give the problem a new name, in this case we call:\n",
    "\n",
    "\\begin{equation}\n",
    "\\zeta_l(t) = \\omega_l \\int_{-\\infty}^t \\epsilon(t') \\cdot \\exp(-\\omega_l (t-t')) dt'\\notag\n",
    "\\end{equation}\n",
    "\n",
    "**memory variables** for the L Maxwell bodies $l=1,...,L$, leading to the stress-strain relation\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma(t) = \\mu_u \\cdot \\epsilon(t) - \\sum_{l=1}^L \\delta\\mu a_l \\zeta_l(t),\\notag\n",
    "\\end{equation}\n",
    "\n",
    "where we do not have to apply the time convolution anymore. However, we need additional partial differential equations in order to estimate the memory variables $\\zeta_l$. So in step 2 we take the first time derivative of the memory variables:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{\\zeta}_l(t) &= \\omega_l \\frac{\\partial}{\\partial t}\\int_{-\\infty}^t \\epsilon(t') \\cdot \\exp(-\\omega_l (t-t')) dt'\\notag\\\\\n",
    "&= \\omega_l \\biggl\\{-\\omega_l\\int_{-\\infty}^t \\epsilon(t') \\cdot \\exp(-\\omega_l (t-t')) dt' + \\epsilon(t)\\biggr\\}\\notag\\\\\n",
    "&= \\omega_l \\biggl\\{-\\zeta_l(t) + \\epsilon(t)\\biggr\\}\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "With this trick, we could get rid of the time convolution, but at the price of L additional partial differential equations to estimate the memory variables, which have to be solved in the FD code:\n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{\\zeta}_l(t) = \\omega_l \\epsilon(t) - \\omega_l \\zeta_l(t)\\notag\n",
    "\\end{equation}\n",
    "\n",
    "To get the partial differential equations in the **stress-velocity** formulation, we only have to take the first time derivative, yielding:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{\\sigma}(t) &= \\mu_u \\cdot \\dot{\\epsilon}(t) - \\sum_{l=1}^L \\delta\\mu a_l \\xi_l(t),\n",
    "\\notag\\\\\n",
    "\\dot{\\xi}_l(t) &= \\omega_l \\dot{\\epsilon}(t) - \\omega_l \\xi_l(t)\\notag\n",
    "\\end{align}\n",
    "\n",
    "with $\\xi_l = \\dot{\\zeta_l}$. To simplify the parametrization, we introduce the dimensionless **anelastic coefficients** \n",
    "\n",
    "\\begin{equation}\n",
    "Y_l = a_l \\frac{\\delta \\mu}{\\mu_u}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "leading to \n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{\\sigma}(t) = \\mu_u \\cdot \\dot{\\epsilon}(t) - \\sum_{l=1}^L \\mu_u Y_l \\xi_l(t).\n",
    "\\notag\n",
    "\\end{equation}\n",
    "\n",
    "## Equations of motion for the viscoelastic 1D SH-problem\n",
    "\n",
    "Finally, we combine the momentum equation with the stress-strain relation to the Equations of motion for the viscoelastic 1D SH-problem. After replacing the simplified notation from the last lesson:\n",
    "\n",
    "\\begin{align}\n",
    "\\sigma &\\rightarrow \\sigma_{yx}\\notag\\\\\n",
    "\\epsilon &\\rightarrow \\epsilon_{yx}\\notag\\\\\n",
    "\\mu_u &\\rightarrow 2\\mu_u\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "we get the conservation of momentum:\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho \\dot{v}_y = \\frac{\\partial \\sigma_{yx}}{\\partial x} + f_y\\notag\n",
    "\\end{equation}\n",
    "\n",
    "and the viscoelastic stress-strain relation:\n",
    "\n",
    "\\begin{align}\n",
    "\\dot{\\sigma_{yx}}(t) &= 2 \\mu_u \\dot{\\epsilon_{yx}}(t) - \\sum_{l=1}^L \\mu_u Y_l \\xi_l(t)\\notag\\\\\n",
    "\\dot{\\xi}_l(t) &= \\omega_l \\dot{\\epsilon_{yx}}(t) - \\omega_l \\xi_l(t)\\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "with the unrelaxed shear modulus:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_u = \\frac{\\mu_0}{1-\\sum_{l=1}^L Y_l}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "and \n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{\\epsilon}_{yx} = \\frac{1}{2} \\frac{\\partial v_y}{\\partial x}\\notag\n",
    "\\end{equation}\n",
    "\n",
    "With these equations, the 1D viscoelastic SH problem is defined and can be solved by the staggered grid FD approach in the final TEW2 lesson."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "- Estimation of the time-domain relaxation function $\\Psi(t)$ from the frequency-domain complex shear modulus of the GMB by inverse Fourier transform\n",
    "- The time-domain viscoelastic stress-strain relation introduces a time convolution, required at each time-step and spatial subsurface point, which would significantly increase the computational costs of the FD modelling\n",
    "- We can get rid of the time convolution by introducing memory variables, at the price of one additional partial differential equations for each Maxwell body\n",
    "- Finally, we compiled the equations of motion for 1D wave propagation in viscoelastic SH media"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
